# Librarieslibrary(tidyverse)library(gmRi)library(scales)library(gt)library(patchwork)library(gtsummary)library(gt)library(sizeSpectra)library(gganimate)library(glue)library(merTools)conflicted::conflict_prefer("select", "dplyr")conflicted::conflict_prefer("filter", "dplyr")# Processing functionssource(here::here("Code/R/Processing_Functions.R"))# ggplot themetheme_set(theme_gmri(axis.line.y =element_line(),axis.ticks.y =element_line(), rect =element_rect(fill ="white", color =NA),panel.grid.major.y =element_blank(),strip.text.y =element_text(angle =0),axis.text.x =element_text(size =8),axis.text.y =element_text(size =8),legend.position ="bottom"))# labels for factor levelsarea_levels <-c("GoM", "GB", "SNE", "MAB")area_levels_long <-c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")area_levels_all <-c("Northeast Shelf", area_levels)area_levels_long_all <-c("Northeast Shelf", area_levels_long)# table to join for swapping shorthand for long-hand namesarea_df <-data.frame(area =c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"),survey_area =c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),area_titles =c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))# Degree symboldeg_c <-"\u00b0C"
Determining Spectra Minimum Sizes
Size spectra typically focus on the descending slope of the size distribution. This doc covers how that minima is selected and the process for doing so, and then the consequence of shifting the minima of the bounded distribution.
Code
# Data used for Wigley estimateswigley_in <-read_csv(here::here("Data/processed/wigley_species_trawl_data.csv")) %>%left_join(area_df)
The general idea is to bin the individual abundances by log2 size-bins, then determine the peak.
The following functions will aid in the labeling of log2 bins:
Code
# Broad Distribution#log2 bins for easy-of-access#' @title Build Log 2 Bin Structure Dataframe#' #' @description Used to build a dataframe containing equally spaced log2 bins for#' size spectra analysis. Contains details on the left and right limits, midpoint, bin width, #' and a text label for the bins. log2bin number ascends with increasing size for eeasy plotting.#'#' @param log2_min #' @param log2_max #' @param log2_increment #'#' @return#' @export#'#' @examplesdefine_log2_bins <-function(log2_min =0, log2_max =13, log2_increment =1){# How many bins n_bins <-length(seq(log2_max, log2_min + log2_increment, by =-log2_increment))# Build Equally spaced log2 bin df log2_bin_structure <-data.frame("log2_bins"=as.character(seq(n_bins, 1, by =-1)),"left_lim"=seq(log2_max - log2_increment, log2_min, by =-log2_increment),"right_lim"=seq(log2_max, log2_min + log2_increment, by =-log2_increment)) %>%mutate(bin_label =str_c(round(2^left_lim, 3), " - ", round(2^right_lim, 3), "g"),bin_width =2^right_lim -2^left_lim,bin_midpoint = (2^right_lim +2^left_lim) /2) %>%arrange(left_lim)return(log2_bin_structure)}#' @title Assign Manual log2 Bodymass Bins - By weight#'#' @description Manually assign log2 bins based on individual length-weight bodymass #' in increments of 1 on the log2 scale. Returns data with bins assigned based on individual#' length-weight biomass#' #' Uses maximum weight, and its corresponding bin as the limit.#'#' @param wmin_grams Catch data prepared for mle calculation, use prep_wmin_wmax#' @param log2_increment Equally spaced increments to use for log 2 bin sizes. Default = 0.5.#'#' @return#' @export#'#' @examplesassign_log2_bins <-function(wmin_grams, log2_increment =1){#### 1. Set up bodymass bins# filter missing weights size_data <- wmin_grams %>%filter(wmin_g >0,is.na(wmin_g) ==FALSE, wmax_g >0,is.na(wmax_g) ==FALSE)# Get bodymass on log2() scale size_data$log2_weight <-log2(size_data$ind_weight_g)# Set up the bins - Pull min and max weights from data available#min_bin <- floor(min(size_data$log2_weight)) min_bin <-0 max_bin <-ceiling(max(size_data$log2_weight))# Build a bin key, could be used to clean up the incremental assignment or for apply style functions log2_bin_structure <-define_log2_bins(log2_min = min_bin, log2_max = max_bin, log2_increment = log2_increment)# Loop through bins to assign the bin details to the data log2_assigned <- log2_bin_structure %>%split(.$log2_bins) %>%map_dfr(function(log2_bin){# limits and labels l_lim <- log2_bin$left_lim r_lim <- log2_bin$right_lim bin_num <-as.character(log2_bin$log2_bin)# assign the label to the appropriate bodymasses size_data %>%mutate(log2_bins =ifelse( between(log2_weight, l_lim, r_lim), bin_num, NA),log2_bins =as.character(log2_bins)) %>%drop_na(log2_bins) })# Join in the size bins log2_assigned <-left_join(log2_assigned, log2_bin_structure, by ="log2_bins")# return the data with the bins assignedreturn(log2_assigned)}#' @title Calculate Normalized and De-Normalized Abundances#'#' @description For binned size spectra estimation we use the stratified abundance divided by the#' bin widths (normalized size spectra). Another way to present the data is to de-normalize, which #' takes those values and multiplies them by the mid-point of the log-scale bins.#' #' min/max & bin_increments are used to enforce the presence of a size bin in the event that #' there is no abundance. This is done for comparing across different groups/areas that should #' conceivably have the same size range sampled.#'#' @param log2_assigned size data containing the bin assignments to use#' @param min_log2_bin Minimum 2^x value for the size spectra being measured (>=)#' @param max_log2_bin Maximum 2^x value for the size spectra being measured (<)#' @param bin_increment The bin-width on log scale that separates each bin#' @param ... Additional grouping factors with which to aggregate on besides the size bins themselves#'#' @return#' @export#'#' @examplesaggregate_log2_bins <-function( log2_assigned, min_log2_bin =0, max_log2_bin =13, bin_increment =1, ...){# Full Possible Bin Structure# Fills in any gaps log2_bin_structure <-define_log2_bins(log2_min = min_log2_bin, log2_max = max_log2_bin, log2_increment = bin_increment)# Capture all the group levels with a cheeky expand()if(missing(...) ==FALSE){ log2_bin_structure <- log2_bin_structure %>% tidyr::expand(left_lim, distinct(log2_assigned, ...)) %>%full_join(log2_bin_structure) }# Get bin breaks log2_breaks <-sort(unique(c(log2_bin_structure$left_lim, log2_bin_structure$right_lim)))# Get Totals for bodymass and abundances log2_aggregates <- log2_assigned %>%group_by(left_lim, ...) %>%summarise(abundance =sum(numlen_adj, na.rm = T),weight_g =sum(wmin_g, na.rm = T),.groups ="drop")# Join back in what the limits and labels are# The defined bins and their labels enforce the size limits log2_prepped <-left_join(x = log2_bin_structure, y = log2_aggregates)#### Fill Gaps with Zero's?? # This ensures that any size bin that is intended to be in use is actually used log2_prepped <- log2_prepped %>%mutate(abundance =ifelse(is.na(abundance), 1, abundance),weight_g =ifelse(is.na(weight_g), 1, weight_g))#### normalize abundances using the bin widths log2_prepped <- log2_prepped %>%mutate(normalized_abund = abundance / bin_width,normalized_biom = weight_g / bin_width,# # Remove Bins Where Normalized Biomass < 0? No!# normalized_abund = ifelse(normalized_abund < 2^0, NA, normalized_abund),# norm_strat_abund = ifelse(norm_strat_abund < 2^0, NA, norm_strat_abund) )# Add de-normalized abundances (abundance * bin midpoint) log2_prepped <- log2_prepped %>%mutate(denorm_abund = normalized_abund * bin_midpoint,denorm_biom = normalized_biom * bin_midpoint)# Return the aggregationsreturn(log2_prepped)}
Using those functions we can assign the bins based on the body-size column in the dataframe we would feed to our spectra slope workflow.
I can’t auto-locate the tallest bin within all of that geom_histogram plot code, so it will be easier to do it outside and animate them using geom_col.
Here are all of them as an animation:
Code
# Get the plotting summarieswig_l2_aggs <- wig_l2 %>%unite("groupvar", area, season, est_year, sep ="XX") %>%split(.$groupvar) %>%map_dfr(function(x){# Get aggregates l2_aggs <-aggregate_log2_bins(x, bin_increment =1)# Locate Tallest tallest <- l2_aggs %>%arrange(desc(normalized_abund)) %>%slice(1) %>%pull(left_lim)# Add info back to aggregates l2_aggs <- l2_aggs %>%mutate(is_tallest =if_else(left_lim == tallest, T, F))# Returnreturn(l2_aggs)}, .id ="groupvar") %>%separate(groupvar, into =c("area", "season", "est_year"), sep ="XX") %>%mutate(area =factor(area, levels = area_levels_long_all),season =factor(season, levels =c("Spring", "Fall")))# # ggplot with animation# animated_plot <- ggplot(wig_l2_aggs) +# geom_col(# aes(# x = left_lim,# y = normalized_abund,# fill = season,# color = is_tallest),# alpha = 0.6,# linewidth = 1.5) +# scale_y_continuous(# trans = transform_log(base = 2),# labels = label_log(base = 2),# breaks = 2^seq(-10,16, 2)) +# scale_x_continuous(# labels = label_math(expr = 2^.x),# breaks = c(0:15)) +# scale_fill_gmri() +# scale_color_manual(values = c("transparent", "black")) +# facet_grid(area ~ season, scale = "free") +# labs(# y = "Normalized Abundance",# x = "Bodyweight (g)",# title = 'Year: {closest_state}', # Dynamic title for year# fill = "Season") +# transition_states(# est_year,# transition_length = 1,# state_length = 5) +# ease_aes('linear') # Smooth transition# # # # # Animate and save# animate(animated_plot, nframes = 300, fps = 10, width = 1000, height = 800)# anim_save(here::here("faceted_tallest_bin_spectra_wigley.gif"), animated_plot)# Show the GIFknitr::include_graphics(here::here("faceted_tallest_bin_spectra_wigley.gif"))
Bodymass spectra abundance peaks
Looking at Individual Years w/ Observable
I like having the ability to look at them one-by-one, but don’t want to make 800 plots or tabs. I’m using this exercise as an excuse to use observable for interactive selections.
Code
# Reformat years as a date - messes with slider# year_avgs <- mutate(year_avgs, est_year = as.Date(str_c(est_year, "-06-01")))# Define data to use for jsojs_define(data = wig_l2_aggs %>%arrange(est_year))# data = filter(wig_l2_aggs, normalized_abund > 1) %>% # arrange(est_year))
Code
viewof yearSlider = Inputs.range( [1970,2019],// If not pulling from the data {label:"Year:",value:1970,// If not pulling from the datastep:1} // Step size for numeric slider)// Create selections for season and area//viewof seasonSelect = Inputs.select(["Spring", "Fall"], { label: "Season" })viewof areaSelect = Inputs.select(["Northeast Shelf","Gulf of Maine","Georges Bank","Southern New England","Mid-Atlantic Bight"], { label:"Area" })// // Filtering Function// filteredData = transpose(data).filter(function(d) {// return yearSlider == d.est_year && d.area === areaSelect;// })filteredSpring =transpose(data).filter(function(d) {return yearSlider == d.est_year&& d.area=== areaSelect && d.season==="Spring";})filteredFall =transpose(data).filter(function(d) {return yearSlider == d.est_year&& d.area=== areaSelect && d.season==="Fall";})
Code
// Plot with FacetsPlot.plot({facet: {data: filteredSpring,x: d => d.area,// Columns: areay: d => d.season// Rows: season//y: { field: "season", domain: ["Spring", "Fall"] } // Custom row order - doesn't work },marks: [ Plot.barY(filteredSpring, { x:"left_lim",y:"normalized_abund",fill:"#00608a",//stroke: "season",stroke: d => d.is_tallest?"black":"none",// Change stroke color for highlighted barsstrokeWidth: d => d.is_tallest?4:0// Thicker stroke for highlighted bars }) ],x: {label:"Body Mass (g)",tickFormat: d =>`2^${d}`// Format labels as powers of 2 },y: {label:"Abundance",type:"symlog",base:2 },style: {facetPadding:5 }})
Code
// Plot with FacetsPlot.plot({facet: {data: filteredFall,x: d => d.area,// Columns: areay: d => d.season// Rows: season//y: { field: "season", domain: ["Spring", "Fall"] } // Custom row order - doesn't work },marks: [ Plot.barY(filteredFall, { x:"left_lim",y:"normalized_abund",fill:"#ea4f12",//stroke: "season",stroke: d => d.is_tallest?"black":"none",// Change stroke color for highlighted barsstrokeWidth: d => d.is_tallest?4:0// Thicker stroke for highlighted bars }) ],x: {label:"Body Mass (g)",tickFormat: d =>`2^${d}`// Format labels as powers of 2 },y: {label:"Abundance",type:"symlog",base:2 },style: {facetPadding:5 }})
Peak Location Frequencies
Without over-doing the peak locations for multi-modal situations, here are the peak locations and their frequencies for each group:
From this dataframe we can strip out the key group information and the minimum size as a lookup table. This can be joined and used as a filter as a preparatory step before fitting the group spectra.
Code
# # Make the key# min_size_key <- wig_l2_aggs %>%# filter(is_tallest) %>%# distinct(area, season, est_year, left_lim, is_tallest) %>%# mutate(min_cutoff_g = 2^left_lim) %>%# select(-left_lim) %>%# mutate(est_year = as.numeric(est_year))# # # # # Join the data back to the original# wigley_in_new <- wigley_in %>%# mutate(area = "Northeast Shelf") %>%# bind_rows(wigley_in) %>%# left_join(min_size_key, join_by("est_year", "area", "season")) %>%# filter(wmin_g >= min_cutoff_g)# # # # # Re-Run the MLE method to get new data# # # and again for 16g to 50kg# peak_chase_results <- group_binspecies_spectra(# ss_input = wigley_in_new,# grouping_vars = c("est_year", "season", "survey_area"),# abundance_vals = "numlen_adj",# length_vals = "length_cm",# use_weight = TRUE,# isd_xmin = NULL,# global_min = FALSE,# isd_xmax = NULL,# global_max = FALSE,# bin_width = 1,# vdiff = 2) %>% # mutate(est_year = as.numeric(est_year)) %>% # left_join(area_df)# # # shelf_peak_results <- group_binspecies_spectra(# ss_input = filter(wigley_in_new, area == "Northeast Shelf"),# grouping_vars = c("est_year", "season"),# abundance_vals = "numlen_adj",# length_vals = "length_cm",# use_weight = TRUE,# isd_xmin = NULL,# global_min = FALSE,# isd_xmax = NULL,# global_max = FALSE,# bin_width = 1,# vdiff = 2) %>% # mutate(# est_year = as.numeric(est_year),# area = "Northeast Shelf") %>% # left_join(area_df)# # # # Join those together for one file# moving_peak_spectra <- bind_rows(peak_chase_results, shelf_peak_results)# # Save them# write_csv(min_size_key, here::here("Data/processed/wigley_species_l2peaks_key.csv"))# write_csv(moving_peak_spectra, here::here("Data/processed/wigley_species_l2peaks_bmspectra.csv"))
# Function to get the Predictions# Flag effect fits that are non-significant ###get_preds_and_trendsignif <-function(mod_x){ modx_preds <-as.data.frame(# Model predictionsggpredict( mod_x, terms =~ est_year * area * season) ) %>%mutate(area =factor(group, levels = area_levels_long_all),season =factor(facet, levels =c("Spring", "Fall")))# Just survey area and year modx_emtrend <-emtrends(object = mod_x,specs =~ area*season,var ="est_year") %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T))# Preds with signif modx_preds %>%left_join(select(modx_emtrend, area, season, non_zero))}
Code
# Model linear trends in timelibrary(nlme)library(ggeffects)library(emmeans)library(tidyquant)moving_peak_spectra <- moving_peak_spectra %>%mutate(area =factor(area, levels = area_levels_long_all),est_year =as.numeric(est_year))spectra_trend_mod <-gls( b ~scale(est_year) * area * season,correlation =corAR1(form =~ est_year | area/season),data = moving_peak_spectra)#check_model(spectra_trend_mod)# plot(check_collinearity(spectra_trend_mod))# Get the predictions and flag whether they are significantbmspectra_trend_predictions <-get_preds_and_trendsignif(spectra_trend_mod) %>%mutate(metric ="bodymass_spectra",area =factor(area, levels = area_levels_long_all))# Make the plotbmspectra_trend_p <-ggplot() +geom_ribbon(data =filter(bmspectra_trend_predictions, non_zero ==TRUE),aes(x, ymin = conf.low, ymax = conf.high, fill = season),linewidth =1, alpha =0.35) +geom_point(data = moving_peak_spectra,aes(est_year, b, color = season),size =0.8, alpha =0.5) +geom_ma(data = moving_peak_spectra,aes(est_year, b, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA, key_glyph ="timeseries" ) +geom_line(data =filter(bmspectra_trend_predictions, non_zero ==TRUE),aes(x, predicted, color = season, linetype ="Significant Trend"),linewidth =1, alpha =1) +scale_fill_gmri() +scale_color_gmri() +facet_grid(area ~ ., scales ="free") +labs(title ="Trends in Spectra w/ Shifting wmin",x ="Year",y ="Exponent of Mass Size Spectra",linetype ="Trend",color ="Season",fill ="Season")# Show plotbmspectra_trend_p
Things to check in the morning: - autocorrelative structure - are seasons correlated now - what is the relationship of steepness to minimum weght cutoff?
Autocorrelation Check
These are annual chek-ins of the multi-species, multi-cohort community size structure in an area. There should be some memory of the previous six months and/or the previous year.
Code
# Pull out acf stuffspectra_acf <- moving_peak_spectra %>%unite("regseas", area, season, sep ="XX") %>%split(.$regseas) %>%map_dfr(function(x){# drop NA x <-arrange(x, est_year) %>%drop_na(b)# Do ACF x_acf <-acf(x$b, plot = F, lag.max =10)# Get the signif:# https://www.squaregoldfish.co.uk/programming/r_acf_significance.md/# Not 100% sure n is the same for ccf as it is for acf, but... significance_level <-qnorm((1+0.95)/2)/sqrt(sum(!is.na(x$b)))data.frame("acf"= x_acf$acf,"lag"= x_acf$lag,"sigpos"= significance_level,"signeg"= significance_level*-1 ) }, .id ="regseas") %>%separate("regseas", into =c("area", "season"), sep ="XX") %>%mutate(# Flag when it crosses thresholdsig_flag =ifelse(acf < signeg | acf > sigpos, T, F),# Set Factor Levelsarea =factor(area, levels = area_levels_long_all),season =factor(season, levels =c("Spring", "Fall"))) # Plot thingsggplot(spectra_acf, aes(lag, acf)) +geom_hline(yintercept =0, color ="gray25", linewidth =1) +geom_vline(xintercept =0, color ="gray25", linewidth =1) +geom_ribbon(aes(ymin = signeg, ymax = sigpos), alpha =0.2, linetype =2, fill ="lightgray",color ="gray65") +geom_col(alpha =0.65) +geom_col(data =filter(spectra_acf, sig_flag), color ="black") +scale_x_continuous(limits =c(-1,10),breaks =c(1:9),expand =expansion(add =c(0,0))) +facet_grid(area~season) +labs(title ="Spectra Slope ACF",x ="Lag (years)", y ="ACF", )
Based on the auto-corrrelation function results we have a similar outlook on the annual autocorrelations.
Spring and fall communities from the same area remain uncorrelated. This is surprising to me as the noise coming from the presence/absence of cohorts of smaller individuals should be reduced.
Code
moving_peak_spectra %>%ggplot(aes(xmin_actual, b)) +geom_point() +geom_smooth(method ="lm") + ggpmisc::stat_poly_eq(ggpmisc::use_label(c("P", "R2")), formula = y ~ x) +scale_y_continuous(expand =expansion(mult =c(0, 0.5))) +scale_x_continuous(trans =transform_log(base =2), labels =label_log(base =2)) +facet_wrap(~area, ncol =2, scales ="free") +labs(title ="Bias from Shifting the Distribution's wmin")
Affinities to Large External Factors
Here I could explore whether the relationships to bottom temperature or landings, but I remain unconvinced this is the right path forward.
# Just drop the random effect part for nowfull_model_ar <- nlme::gls(# Modelmodel = b ~ area * season *scale(bot_temp) + area *scale(log(total_weight_lb)), # Datadata = peak_chase_model_df %>%drop_na(total_weight_lb) %>% dplyr::filter(area !="Northeast Shelf"),# Autocorrelationcorrelation =corAR1(form =~ est_year | area/season))# Check the modellibrary(performance)# check_model(full_model_ar)# # Collinearity without interactions# plot(# check_collinearity(# nlme::gls(# b ~ area + season + scale(bot_temp) + scale(log_land),# data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>%# filter(area != "Northeast Shelf") %>% # mutate(log_land = log(total_weight_lb)),# correlation = corAR1(form = ~ est_year | survey_area/season)# )# )# )# # confidence intervals on phi# intervals(full_model_ar)$corStruct
Code
# Clean up the plot:# Plot the predictions over datatemp_preds <-as.data.frame(ggpredict( full_model_ar, terms =~ bot_temp * area * season)) %>%mutate(area =factor(group, levels = area_levels_long_all),season =factor(facet, levels =c("Spring", "Fall")))#### Trend Posthoc ##### Flag effect fits that are non-significant ###temp_emtrend <-emtrends(object = full_model_ar,specs =~ area * season,mode ="appx-satterthwaite",var ="bot_temp") %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T))# Limit temperature plotting rangeactual_range <- peak_chase_model_df %>%filter(area !="Northeast Shelf") %>%group_by(season, area) %>%summarise(min_temp =min(bot_temp)-2,max_temp =max(bot_temp)+2,.groups ="drop")# Plot over observed datalocal_btemp <- temp_preds %>%left_join(actual_range) %>%filter((x < min_temp) == F, (x > max_temp) == F) %>%left_join(select(temp_emtrend, area, season, non_zero)) %>%filter(non_zero) %>%ggplot() +geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha =0.1) +geom_line(aes(x, predicted, color = season), linewidth =1) +geom_point(data =filter(peak_chase_model_df, area!="Northeast Shelf"),aes(bot_temp, b, color = season),alpha =0.4,size =1) +facet_grid(area~., scales ="free", labeller =label_wrap_gen(width =8)) +scale_color_gmri() +scale_x_continuous(labels =label_number(suffix = deg_c)) +labs(y ="Mass distribution Exponent",x ="Local Seasonal Bottom Temperature Anomaly")local_btemp
Code
# Clean up the plot:# Plot the predictions over dataf_preds <-as.data.frame(ggpredict( full_model_ar, terms =~ total_weight_lb * area * season)) %>%mutate(area =factor(group, levels = area_levels_long_all),season =factor(facet, levels =c("Spring", "Fall")))#### Trend Posthoc ##### Flag effect fits that are non-significant ###f_emtrend <-emtrends(object = full_model_ar,specs =~ area,var ="total_weight_lb",mode ="appx-satterthwaite" ) %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T))# Limit temperature plotting rangeactual_range <- peak_chase_model_df %>%group_by(season, area) %>%summarise(min_f =min(total_weight_lb)-2,max_f =max(total_weight_lb)+2,.groups ="drop")# Plot over observed datalandings_ar_plot <- f_preds %>%left_join(actual_range) %>%filter((x < min_f) == F, (x > max_f) == F) %>%left_join(select(f_emtrend, area,non_zero)) %>%filter(non_zero) %>%ggplot() +geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha =0.1) +geom_line(aes(x, predicted, color = season), linewidth =1) +geom_point(data =filter(peak_chase_model_df, area !="Northeast Shelf"),aes(total_weight_lb, b, color = season),alpha =0.4,size =1) +facet_grid(area~., scales ="free", labeller =label_wrap_gen(width =8)) +scale_color_gmri() +scale_x_log10(labels =label_log(base =10), limits =10^c(0,10)) +labs(y ="Exponent of ISD",title ="Total Landings (lb.)",x ="Local Seasonal Bottom Temperature Anomaly")landings_ar_plot
Bottom Temperature CCF
Shelf-Wide:
Regionally:
Code
# Function to grab the correlation data and lag dataget_ccf_vector <-function(x,y, lagmax =10){# Run the ccf ccf_data <-ccf(x, y, plot= F , lag.max = lagmax)# Get the signif:# https://www.squaregoldfish.co.uk/programming/r_acf_significance.md/# Not 100% sure n is the same for ccf as it is for acf, but... significance_level <-qnorm((1+0.95)/2)/sqrt(sum(!is.na(x)))data.frame("acf"= ccf_data$acf,"lag"= ccf_data$lag,"sigpos"= significance_level,"signeg"= significance_level*-1 )}# Run Local CCF on Bottom Temperaturebtemp_ccf <- peak_chase_model_df %>%unite(col ="regseas", area, season, sep ="XX") %>%split(.$regseas) %>%map_dfr(~get_ccf_vector(x = .x$bot_temp, y = .x$b), .id ="regseas") %>%separate(col ="regseas", into =c("area", "season"), sep ="XX") %>%mutate(# Flag when it crosses thresholdsig_flag =ifelse(acf < signeg | acf > sigpos, T, F),# Set Factor Levelsarea =factor(area, levels = area_levels_long_all),season =factor(season, levels =c("Spring", "Fall")))# Plot itggplot(btemp_ccf, aes(lag, acf)) +geom_hline(yintercept =0, color ="gray25", linewidth =1) +geom_vline(xintercept =0, color ="gray25", linewidth =1) +geom_ribbon(aes(ymin = signeg, ymax = sigpos), alpha =0.2, linetype =2, fill ="lightgray",color ="gray65") +geom_col(alpha =0.65) +geom_col(data =filter(btemp_ccf, sig_flag), color ="black") +scale_x_continuous(expand =expansion(add =c(0,0))) +facet_grid(area~season) +labs(title ="Bottom Temperature CCF",x ="Bottom Temperature Lag (years)", y ="ACF", )
Landings CCF
Its plausible that landings have a lagged effect so here are the CCFs for that.
Code
# Run Local CCF on Bottom Temperaturelandings_ccf <- peak_chase_model_df %>%drop_na(total_weight_lb) %>%unite(col ="regseas", area, season, sep ="XX") %>%split(.$regseas) %>%map_dfr(~get_ccf_vector(x = .x$total_weight_lb, y = .x$b), .id ="regseas") %>%separate(col ="regseas", into =c("area", "season"), sep ="XX") %>%mutate(# Flag when it crosses thresholdsig_flag =ifelse(acf < signeg | acf > sigpos, T, F),# Set Factor Levelsarea =factor(area, levels = area_levels_long_all),season =factor(season, levels =c("Spring", "Fall")))ggplot(landings_ccf, aes(lag, acf)) +geom_hline(yintercept =0, color ="gray25", linewidth =1) +geom_vline(xintercept =0, color ="gray25", linewidth =1) +geom_ribbon(aes(ymin = signeg, ymax = sigpos), alpha =0.2, linetype =2, fill ="lightgray",color ="gray65") +geom_col(alpha =0.65) +geom_col(data =filter(landings_ccf, sig_flag), color ="black") +scale_x_continuous(expand =expansion(add =c(0,0))) +facet_grid(area~season) +labs(title ="Total Landings CCF",x ="Total Landings Lag (years)", y ="ACF", )